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The presence of fluctuating local relaxation times, Tp{t) has been used for some time as a con- 
ceptual tool to describe dynamical heterogeneities in glass-forming systems. However, until now no 
general method is known to extract the full space and time dependent Tp(t) from experimental or nu- 
^— ■>) , merical data. Here we report on a new method for determining the local phase fleld, 4>r{t) = /' 

^Sj ■ from snapshots {f{ti)}i=\,.,M of the positions of the particles in a system, and we apply it to extract 

I ^ I <l>r{l) and rff(t) from numerical simulations. By studying how the phase fleld depends on the number 

■ of snapshots, we find that it is a well defined quantity. By studying fiuctuations of the phase field, 

■ we find that they describe heterogeneities well at long distance scales. 

I • I. INTRODUCTION 

' As glass-forming liquids enter the supercooled state, relaxation processes gradually but dramatically slow down, 

c/j . and at the same time they become non-exponential [l|-l3| . One possible way of understanding this behavior would be 

■ to think of the relaxation function F{t) as the sum of a large number of exponentially decaying contributions, each 
^ ; , with different relaxation times; with each contribution corresponding to a particular region in the system 0, This 

■ is one aspect of the picture of dynamical heterogeneity , in which mesoscopic regions relax differently from each other 
and from the bulk [ll4l2l| . In this framework, it is assumed that the differences between regions are not frozen in but 
dynamical, i.e. as the system evolves, the local relaxation times (t) evolve, so that "fast" regions become "slow" , 

' O ' and viceversa. This same picture has also been invoked to explain other phenomena occuring in the same regime, 
C . such as the breakdown of the Stokes-Einstein relation between the diffusion coefficient and the shear viscosity. Direct 
^ ■ evidence for dynamical heterogeneity has been found in particle tracking experiments in glassy colloidal systems 
and granular systems @ , and in numerical simulations [l^ [ll| . 

The origin of dynamical heterogeneity is still uncertain, although some explanations have been proposed [T^ . One 
of the proposed explanations is that dynamical heterogeneity emerges as a consequence of dynamical constraints, 
which produce non trivial structure in the space of (space-time) particle trajectories [H, Ell- Another proposal is 
^ _ that as a liquid is supercooled it eventually undergoes a Random First Order Transition in which the liquid freezes 
. into a mosaic of aperiodic crystals [l^ [l^, which give rise to the heterogeneity that is observed. A third approach 
postulates that dynamical heterogeneity emerges from the presence of soft (Goldstonc) modes associated with a broken 
continuous symmetry under time reparametrizations t — ^ 4>{t) p7l - l27j . 

A quantitative description of dynamical heterogeneity in terms of the presence of locally fluctuating relaxation 
^vq _ times Tp{t) has in principle some strong advantages. One of them is its simplicity and intuitive appeal. Another one is 
T-H ■ that the basic quantities that appear in this description are intrinsically instantaneous, as opposed to other common 
\ descriptions for which the basic quantities describe changes in the system over a finite time interval, which complicates 
their interpretation. Despite those advantages, however, such a description has proved elusive. For example, even the 
question of experimentally determining the relationship between the lifetime of regions of heterogeneous dynamics 
JH ■ and the bulk a-relaxation time Tq, has proved controversial, with some results indicating that Tc^/Ta ^ 1, and other 
results indicating that T^-^jTa 3> 1 

In the present work, our goal is twofold. On one hand, we establish a connection between the presence of time 
reparametrization fluctuations and the presence of fluctuating local relaxation times Tp{t). On the other hand, we 
introduce a new method to extract the actual values of those relaxation times Tp(t) from experimental or numerical 
data, and test this method on data from numerical simulations. 

Our paper is organized as follows. In Sec.|TT]we discuss the connection between the presence of time reparametriza- 
tion fluctuations and a description of the system in terms of fluctuating local relaxation times. In Sec. IIIII we explain 
our proposed method to extract local relaxation times from numerical or experimental data. In Sec. IIVI we discuss 
the details of the numerical simulations that we have used to test the method, and we present the results of those 
tests. Finally, in Sec. [Vj we summarize our results and present some remarks and ideas for further work. 
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II. CONNECTION BETWEEN FLUCTUATING LOCAL RELAXATION TIMES AND TIME 

REPARAMETRIZATION FLUCTUATIONS 



At temperatures close to the calorimetric glass transition temperature Tg, the dynamics of glass forming liquids 
is dramatically different from the high temperature dynamics 0, [28| . At higher temperatures the relaxation is 
exponential, a direct consequence of Brownian-like motion. By contrast, at lower temperatures the relaxation is 
exponential only at the early times (with fast dynamics) and is followed by a plateau at intermediate times. This 
plateau is attributed to the caging effect that results from having each particle temporarily trapped by its neighbors. 
As the cages start to break and the particles move distances of the order of the intcr-particle distances, the relaxation 
function starts to decay. Unlike what happens at high temperatures, the final decay is usually not exponential and in 
many cases it is described well by the Kohlrausch- Williams- Watts (KWW) "stretched exponential" function 



C{t-t')=Aexp{~[{t~t')/Tf}, 



(1) 



where we are using some correlation function C{t — t') between the states of the system at times t and t' as an example 
of a relaxation function. Here A and /3 are constants and (3 is sometimes called the "stretching exponent" or the 
" Kohlrausch exponent" . Two extreme scenarios were proposed to explain the non-exponential behavior [2I, '28*1 . In one 
of them the liquid is heterogeneous, the relaxation is exponential in each small region, there is a spatial distribution 
of relaxation times P{t), and the relaxation function is given by 



C{t - t') = 



f°° f t-i'\ f (fir 
J drP(r)exp f — j = J —Cp{t,t') with Cf{t,t') = Co cxp 



t-t' 



(2) 



where Cf(i, t') is the local two-time correlation in a small region around point r and is the corresponding relaxation 
time. In the other extreme scenario, the dynamics of the liquid is homogeneous, the distribution of the relaxation times 
is a delta function, and the relaxation in all regions is described by a unique function g{x), which is non-exponential: 



C{t - 1') 



drSir - To)g 



t-t' 



Cr{t,t') with. Cp{t,t') =^ g 



t-t' 
To 



(3) 



Eqs. ^ and ^ can be generalized to allow for the simultaneous presence of both heterogeneous relaxation times and 
non-exponential local relaxation functions [29j : 



C{t-t')= / dTPg{T)g 



t-t' 



= / 



with Cf{t,t') = g 



t-t' 



(4) 



where Pgir) is the probability density of relaxation times compatible with the local relaxation function g{x). 

If we allow the time difference t — t' to be long enough, we cannot assume that each local region is characterized by a 
single relaxation time rp, because rp could in principle have fluctuated over this time interval. To take this effect into 
account, let us divide the time interval [t',t) into n — l sub-intervals {[^i-i, ii)}i=2,...,n such that for each sub-interval 
the fluctuations in Tp{t) are negligible. Hence, for long time intervals, the local correlation is given by 



cAt,t')^, 



n 



(5) 



where tn 
integral. 



t and ti = t' . In the limit of infinitesimal time intervals we get ti — ti-i — ?■ dt and the sum becomes an 



Cr{t,t') 



dt" 



Mt") 



(6) 



In this general expression for Cp{t, t'), the relaxation is allowed to be locally non-exponential, and the relaxation time 
is allowed to fluctuate both in space and time. 

Another way of looking at local correlations is to consider first the exponential relaxation case where 



dC{t,t') 
Jt 



-C{t,t') 



C{t,t') = Cocxp 



t-t' 



(7) 



If we allow r to fluctuate in space and time but postulate that the same differential equation still holds, then we have 



dCp{t, t' 
Jt 



1 



Tp{t) 



Cr{t,t') 



Cp{t, t') = Co cxp 



df' 



Mt") 



(8) 
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where Cp{t,t') and Tf{t) are defined as before and we have used the initial condition Co = limj_^t/+ C{t,t'). Eq. ^ 
can be extended to describe the more general case of non-exponential relaxation. For example, we can modify the 
right hand-side of the differential equation in Eq. ([S]) to = —^G[Cp{t,t')]. Here the function G[Cp{t,t')] 

represents the effect of non-exponential relaxation. The solution of the new differential equation is 

Cr-{t,t')^c[Mt)-Mt')], (9) 

where (l)p{t) = J* t^tj, and C{x) is a monotonous decreasing function that satisfies C'{x) = —G[C{x)] and < C{x) < 
1. 

By comparing Eq. ^ with Eq. (|9]) we find that g{x) — C{x) and 

Mt)-Mt')^ (10) 

A completely different way of deriving Eq. 1^ is based on the presence of time reparametrization symmetry, and we 
describe it next. Time reparametrization symmetry refers to an invariance under transformations of the time variable 
t 4){t), where (t){t) is a continuous and monotonically increasing function of t. This symmetry has been shown to 
be present in mean field spin glass models [sol . [3l| a.nd in short range spin glass models [13, [H [HI, HI] . Numerical 
studies in spin glasses [19|] and structural glasses |2"ol. [26| have given strong evidence supporting the presence of this 
symmetry. The symmetry has also been used, from a different point of view, to determine properties of symmetric 
observables whose full dynamical properties are hard to compute [s^ [33| . The symmetry is spontaneously broken by 
correlations and relaxation functions. We can illustrate the breaking of the symmetry by considering the correlation 
function C{t — t'). If C{t — t') was symmetric under time reparametrizations t — ?► (^(t), then C{t — t') ~ C[(f>{t) — (f>{t')] 
for all mappings t — which is only possible if the correlation is independent of time. But we know that in the 
long time limit the correlation function of a glass forming liquid decays from its plateau value, therefore the symmetry 
is broken by C{t — t'). The broken reparametrization symmetry is expected to give rise to Goldstone modes [17H27| 
associated to smooth variations in space of the reparametrized variable, that is t — >■ <j)r{t)- The spatio-temporal 
fluctuations of imply that different regions of a sample relax differently from each other, that the relaxation 

of the different regions is advanced (or retarded) with respect to each other and the bulk, and that "advanced" and 
"retarded" regions can switch roles in the course of the relaxation. Hence, the spatio-temporal fluctuations induced by 
the broken reparametrization symmetry provide a possible explanation for dynamical heterogeneities in glass forming 
and glassy systems. In addition, the broken reparametrization symmetry leads to recovering Eq. Q if we ignore 
longitudinal fluctuations tlZ-lft] . Longitudinal fluctuations should be suppressed by coarse graining in larger regions 
because at long distances they are less correlated than transverse fluctuations [13] ■ Finally, it should be noticed that 
there is a gauge symmetry under time-independent shifts of the reparametrization variable (t) , that is in Eq. 1^ , 
the local correlation is unchanged by the transformation (j)p{t) — ^ 4'p{t) + pp. Consequently, when we analyze data we 
cannot determine an absolute (prit)- To work with gauge- invariant quantities, we either have to go back to studying 
two-time quantities like 0f7(i) — (j)r{t') or we must work with time derivatives of 4ir{i)- 



III. METHOD TO EXTRACT PHASES AND INSTANTANEOUS RELAXATION RATES FROM 

NUMERICAL OR EXPERIMENTAL DATA 

Our method for the extraction of (j)r{t) from local two-time correlations is inspired by the work in Rcfs. [l^ll^l. The 
authors of Refs. [l^H^] use a functional form 5giobai(2;) = qEA exp (— |a;|'') , where qsA and /3 are fltting parameters, to 
fit global correlations in molecular dynamics simulations of particle systems and polymer systems. They find functions 
(t){t) for the different systems such that the data are described by 

C7(i,t') = .9giobai[<^W "0(i')], (11) 
As discussed in Sec. [Hi local fluctuations can be described in terms of Goldstone modes such that 

Cr{t,t')^g[cpr{t)~^p{t')]. (12) 

In principle, the functions (/global (a^) and g(x), respectively describing global two-time correlations and local two-time 
correlations, could have completely different forms. For example, as discussed in Sec. [Hi some of the initial motivation 
for the picture of dynamical heterogeneity came from the idea that non-exponential relaxation in a macroscopic 
sample is due to the combined effect of local exponential relaxations with different relaxation times. In this picture, 
5giobai(a:) cx exp(— |a;|'') and g{x) oc exp(— |a;|). In Refs. as a simplifying assumption, the two functions were 
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taken to be equal. In the present work, we focus exclusively on the local function g{x). However, we will consider 
different coarse graining sizes when determining Cp{t, t'). Each of these coarse graining sizes will give rise to a different 
form for g{x). In practice, to simplify the determination of g{x), we will restrict it to be a member of a family of 
functional forms g{x]a)^ parametrized by a vector a of p components, and for each coarse graining size both the 
parameter vector a and the fluctuating phases (t>p{t) will be determined by fitting Cp{t, t') according to Eq. (fT2)) . For 
technical reasons, we impose the condition that g be an even function, g{—x) = g{x). 

Let's consider a data set containing "snapshots" of the relevant degrees of freedom of the system, taken at times 
{ti\i=i,...,M- In our case the positions of the particles are recorded at each time step. For a given coarse graining 
region centered at a point r in the sample, our data points will be the M{M — l)/2 two-time local correlations 
Cp(ti,tj) calculated from the recorded data. Considering for the moment a fixed parameter vector a, we have M 
fitting parameters {^F(ti)}i=i,...,Af j thus giving us a fitting problem with ~ M/2 data points per fitting parameter. 
However, since the function g(x) is nonlinear, this becomes a nonlinear fitting problem with a large number M of 
fitting parameters. In order to make this problem manageable, we want to convert it into a linear fitting problem. To 
do this, we expand the rhs of Eq. (fT2|) to linear order in the local fluctuations 

g[Mt) - Mt')] ~ gm) - Ht')] + g'm) - [-^0^(0 - SMt')] (i3) 

where (j>{t) is a global phase, still to be determined, and the local fluctuations of the phase are given by 

6Mt) = Mt) - m- (14) 

At this point, we have simplified the problem by converting it into a linear fitting problem, at the price of introducing 
the extra phase variables {(/'(ti)}i=i,... ,m. 

Before we describe the method for determining the phases, we define three quantities to measure fluctuations in 
the two-time local correlations. The total fluctuations are defined by 

5Cp{t, t') EE a-(i, t') - gm) - 0(i')] (15) 

and contain complete information about the space dependence of the local two time correlations. We decompose the 
total fluctuations as the sum of a transverse component 5Cp{t,t') and a longitudinal component 6Cp{t,t'), 

5Cp{t,t') = SCj{t,t') + 6Ci{t,t'). (16) 

The transverse component is associated to time reparametrization fluctuations, and it is defined by the linear term 
in Eq. (fT3)) . namely 

SC^it, t') ^ g'm) - m] [5Mt) - Hr{t')] . (17) 

To satisfy Eq. the longitudinal component is defined by 

SCi{t,t') = SCp{t,t')~SCj{t,t') (18) 
= Cpit,t') ~ {gm - 0(i')] + g'm) - Ht')] [SMt) - 6Mt')]} ■ (19) 

In what follows we explain how to perform a least squares fit of the local correlation Cf (t, t') by the rhs of Eq. (|13p . 
In the case of numerical data, we normally have not only data for different coarse grained regions inside the sample, 
but also data for different independent runs. If wc have Nn independent simulation runs for a system of volume V, we 
can imagine juxtaposing the configurations for all the runs, thus creating a larger volume V = NfiV. Assuming that 
each coarse graining region Bp centered around point f has a volume Vcg, the total number of non-overlapping regions 
per run is w = [V/Vcg], and the total number of non-overlapping regions including all runs is f2 = Nj^uj — Nii[V/Vcg], 
where in this context [V/Vcg] denotes the integer part of V/Vcg- 

The fit we want to perform corresponds to minimizing the quantity 

E= -^'^e[{S(f>r{ti)}i=i,--- ,M;a,{Cp{tj,ti)}i<i<j<M] (20) 
f 

with respect to all phase fluctuations and with respect to a. Here 

e[{^M^^)}^=l,•••,M;c^,{a-(^J,^^)}l<^<,<Af] ^ry^^Af) J2 [S0{tj,t,)]\ (21) 

l<i<j<A/ 

= V-HM) iCAt,t,)-{g[,pitj)~,pit,);a]+g'm,)~,pit,);a][SMtj) (22) 

l<i<j<A/ 
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is the residual corresponding to the region centered at r, 

,(M) . - M + 1 ^ MiM-3) + 2 ^^3^ 

is the number of degrees of freedom in the fit of the data corresponding to one particular coarse graining region, and ^ - 
denotes a sum over non-overlapping coarse graining regions. From now on we consider only times belonging to the set 
{ti}i=i,--- ,M, and we simplify our notation by indicating times as subindices, i.e. 6(f>if = 5(t>p{ti), SC^^- = SCi'{tj,ti), 
and similarly for all other one- and two-time variables. 

To minimize E, we observe that for a fixed parameter vector a, the determination of the local fluctua- 
tions {6(j)ip} for a specific region centered at f can be performed by separately minimizing the corresponding 
,M] Oi, {Cjip}i<i<j<M) for that particular region. For this reason, instead of minimizing with respect to 
all variables in one step, the problem becomes significantly simplified if we perform an iterated minimization, 

min E — mine0(cJ), (24) 
e^(d) = ^V" min e ({50^?}; a, {Cjif?}i<i<j<Af ) . (25) 

f 

In other words, we first keep a fixed and separately minimize with respect to the phase fluctuations in each coarse 
graining region, and then minimize the result with respect to a. The minimization with respect to the phases, at 
fixed (J, is performed in two steps. In the first step we determine global phases. In the second step, we use the global 
phases from the first step to determine local phases. Next, we describe these steps in detail, but without derivations. 
The details of the derivations can be found in Appendix [Xj 



A. Step One: determining global phases at fixed a 

The global phases {4>i\i=i,...,M define the point around which the Taylor expansion in Eq. ([T^ is performed. The 
only requirement on them is to be close enough to the values of (j)ifr such that the expansion is accurate to first order 
in the fiuctuations. In particular, we only need to determine them to within an error of the order of the fluctuations 
in the 4>ip. There is more than one possible way to satisfy these conditions, and for simplicity we do it by choosing 
them to be the phases that best represent the global correlations for the given value of a. We define the quantity e 

e-({0i, • • • , c^mY a) = r\M) ^ [c,, ~ g'^^\c^, - a)] ' , (26) 

l<i<j<M 

where Cji = C{tj, ti), and g^^\(f)j — (f>i; a) is the first order Taylor expansion of g{(l)j — (f)i\a) with respect to the phase 
difference, taken about 

= [sgxv{t,-U)]g-^[C{t„U)]. (27) 

To minimize e we impose the condition that all derivatives of e with respect to (pk must be zero. Additionally, we fix 
the gauge with the condition 



Thus we obtain the matrix equation 



where 



= ^0,. (28) 



w = A(j) (29) 



_ _ f for fc = 1 

= iE;ii<i>...9'^(<i>..;«) for Ml ^''^ 

_ r 1 for fc = 1 

= I -4.E^iiff"(*.fc;") + .9"(*fc.;«) for fc^l. ^^^^ 

The solution to Eq. (j29p is the set of global phases {(/)i}i=i,...,M which minimizes e for the given d. We use the global 
phases as input in the second step, which we describe next. 
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B. Step Two: determining the local phases at fixed a 



In the second step, we minimize t{{5(j)ip\] a) with respect to the local phase fluctuations {5(j)ift, • • • , 5(j)Mr\i for each 
coarse graining region, while keeping a fixed at the values used in the first step. By Taylor expanding g{(j)ip — 4>jp] a) 
about the global phase differences, 0^ — 4>j, imposing the condition that all derivatives of e{{5(j)ip\] a) with respect to 
the phase fluctuations {54>ip, • • • , 54>Mr} should be zero, and fixing the gauge with the condition 



M 



(32) 



i=i 



we obtain the matrix equation 



where 



wp — A S(f>p, (33) 



_ r for fc = 1 

Wkr^ = I _ ^) j^^. ^ (34) 

_ r 1 for fc = 1 

^ I 4. E-Li - 0, ; a) - - 4; «) for k + 1. ^^^^ 

The solution to this matrix equation is the vector i5^f = ((5(/)if , ■ • • , S4>Mr) containing the local fluctuating phase 
differences that minimize e for the given value of the parameter vector a and for the region Bp. By averaging the 
minimum values of e{{S(j)if}; a) for all regions, we obtain e^{d). 



C. Minimization with respect to a 

The optimum value of a is obtained by numerically minimizing e<;/,(Q?) with respect to it. Once the optimum a has 
been obtained, the solutions of Eq. (|33p for that specific value of a and for each region Br give the optimal choice for 
the fluctuating phases S(f)ip. 



IV. TESTING THE METHOD WITH DATA FROM NUMERICAL SIMULATIONS OF 

GLASS-FORMING MODEL SYSTEMS 



We tested the method on the results of classical Molecular Dynamics simulations of two glass forming systems, each 
one consisting of an 80 : 20 binary mixture of 1000 particles interacting via purely repulsive Weeks-Chandler- Andersen 
(WCA) potentials. The systems were quenched from an initial temperature Ti much higher than the Mode Coupling 
critical temperature Tc- One of the systems was quenched to a temperature T below the mode coupling critical 
temperature Tc such that T/Tc ^ 0.9 and it did not reach equilibrium during the simulation. We call this system the 
aging system. The other system was quenched to a temperature T such that T/Tc 1.1 and it reached equilibrium 
within the duration of the simulation. In our studies we use the data corresponding to times after the system reached 
equilibrium and we call this system the equilibrium system. After the quench, the systems were allowed to evolve 
for times much longer than the typical vibrational periods of the particles. Snapshots of the systems were taken 
periodically during the evolution and the positions of the particles were stored. 

To probe a system, we divide it into coarse graining boxes Bp, centered at the space points f. The local correlation 
for each box is defined by 

Cp{t,t')^^ E cos{q-[nit)-nit')]}- (36) 

Here Nb^, is the number of particles in the coarse graining box Bf at time t' , and \q\ is chosen to be the value 
given by the location of the main peak of the static structure factor S{q). Cp{t, t') measures the extent to which the 
configuration in the coarse graining region has changed between the times t' and t. If the particles initially in the 
region have not moved much then Cp{t,t') is larger than 1/2, and if the particles have moved a significant distance 
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then Cf{t, t') <C 1. In the first case, we say the the region is "slow" , and in the second, that it is "fast" . To implement 
the method, we choose the fit function to be the stretched exponential 

g{x;d) ^ qEAexp{~\x\'^) . (37) 

In this case the vector of parameters is a = {qEAiP)i where qea is the plateau value of the correlation function 
and (3 is the stretching exponent. We begin the analysis by determining the optimal values of a at different coarse 
graining sizes; the results are summarized in Fig. ^ for both systems. Each one of the curves is drawn using two 
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FIG. 1: Variation of the optimal values for the fitting parameters /3 and qsA with coarse graining size, in the equilibrium system 
with M = 94 snapshots and Nr — 10 independent runs (right panel), and in the aging system with M — 100 and A''^ = 10 (left 
panel). The part that is drawn in black for each curve corresponds to the range of coarse graining sizes for which longitudinal 
fluctuations account for more than 40% of the total fluctuations. 



colors: one part is drawn in black and the other in a color other than black. The part that is drawn in black is the 
part where longitudinal fluctuations, i.e. the fit residuals, account for more than 40% of the total variance of the 
fluctuations. We observe that the optimal /3 decreases with increasing coarse graining volume and the optimal qea 
increases with coarse graining volume. Both quantities approach constant values at high coarse graining volumes. 
The trend on qea appears to be an artifact of the fact that fluctuations are stronger for smaller coarse graining 
volumes. Since fluctuations that increase Cp{t, t') above a value given by the Debye- Waller factor are unlikely, there is 
a bias for fluctuations to reduce, rather than increase, the value of Cf{t,t') for times when Cf{t,t') is high. Since the 
largest values that Cp{t,t') takes are associated to the optimal Qea, this bias will push qea down, particularly when 
fluctuations are stronger due to the small size of the coarse graining regions. The trend on /3 is reminiscent of one of 
the initial motivations for the proposal that dynamics must be heterogeneous, i.e. the idea 0, 0] that in systems for 
which the bulk relaxation is non-exponential, local regions have exponential relaxation - /? w 1 - with heterogeneous 
timescales, and this translates into a lower exponent /? for the system as a whole. However, as we will discuss below, 
for smaller coarse graining volumes the majority of the fluctuations are captured by the fit residuals and not by the 
fitting function itself, so it is unclear whether or not this trend in the optimal value of the fit parameter /3 has any 
physical significance. 



A. Fluctuations 



As we mentioned in Sec. [TT] we expect longitudinal fluctuations to be less correlated at long distances, and conse- 
quently to be more strongly suppressed by coarse graining, than transverse fluctuations [s^l- It is not immediately 
obvious how to quantify the strength of the fluctuations, since they are functions of the position and of two times. 
However, in this context it becomes natural to think of them as vectors in an Euclidean vector space, with the inner 
product defined by 

^^^^^ ^ :jM(jr^) ^ Y.(fr<t^^^Mu,t,)), (38) 

^ ' l<j<i<M r 

and the Euclidean norm defined by 



ll/ll^(/l/)' 



(39) 
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Here M is the number of snapshots, {r} are the centers of each one of the to non-overlapping coarse-graining boxes Bp 
in the volume V of the system, and (• • ■ ) denotes an average over thermal fluctuations, which in the case of quantities 
obtained from numerical simulations corresponds to an average over independent simulation runs. The only slightly 
unusual aspect in these definitions is the prefactor 2/ujM{M — 1), which was included so that the square of the norm is 
an "intensive" quantity, which gives the average over positions and over time pairs of the variance of the fluctuations 
for one coarse graining box and one time pair. This makes it easier to compare results for different snapshot numbers 
M and different numbers w of coarse graining boxes. With these definitions, and as a direct consequence of the 
minimization conditions that determine the values of the phases {(f>r{ti)}j we show in Appendix IB] that the transverse 
and longitudinal fluctuation vectors are orthogonal, i.e. 

{6C^\SC^)=0, (40) 

and therefore they satisfy the Pythagorean condition 

\\SC\\^ ^\\SC^\\^ + \\SC^\\\ (41) 

In Fig. ([2]) we show how the strengths of the fluctuations vary with the average number of particles A''b_ per coarse 
graining box. In the figure, the two top panels correspond to the aging system and the two bottom ones correspond 
to the equilibrium system. For each system, the left panel shows the magnitudes ||<5C||, ||^C^||, and ||<5C^|| of the 
total fluctuations, their transverse component and their longitudinal component; and the right panel shows the ratios 
|JC"^||/||5C||, and ||(5C^||/||(5C||. Due to the Pythagorean condition, the sum of the squares of these two ratios is 
unity in all cases. 

We expect that as the coarse graining is increased, all fluctuations will be suppressed by the effect of averaging. 
Indeed, we find that the magnitudes of all three kinds of fluctuations decrease as we increase the coarse graining 
volume, both for the aging system and for the equilibrium system. The magnitude of the longitudinal fluctuations 
decreases at a faster rate than the total fluctuations while the magnitude of the transverse fluctuations decreases at a 
slower rate. Indeed at small coarse graining volumes the ratio | |5C^| |/| |5C| | is smaller than the ratio ||5C^||/||(5C||, 
but the first one increases and the second one decreases as the coarse graining volume increases. For the aging system 
the two ratios cross at Nb^, ~ 30 and for the equilibrium system they cross at A^b_ « 100. This is consistent with 
our expectation that longitudinal fluctuations should be more strongly suppressed by coarse graining than transverse 
fluctuations. We also observe that transverse fluctuations are dominant for a wider range of coarse graining sizes in 
the case of the aging system than in the case of the equilibrium system. This is indeed what should be expected if the 
transverse fluctuations are the Goldstone fluctuations associated with time reparametrization symmetry. Since time 
reparametrization symmetry is a long time asymptotic effect, and in the case of the equilibrium system, which is at 
higher temperature than the aging system, the relaxation time provides a cutoff for long timescales, we expect that 
the transverse fluctuations associated with this symmetry will manifest themselves less strongly in the equilibrium 
system than in the aging one. 



System 




QEA 


/3 


Equilibrium 


216 


0.65 


0.775 


Aging 


125 


0.775 


0.875 



TABLE I: Fitting parameters for coarse graining sizes such that ||(5C"^|p ^ 0.6 ||(5C|p in tire equilibrium system with AI = 94 
and Nr = 10, and in the aging system with AI = 100 and Nr — 10. 

From now on, we consider in detail one coarse graining volume for each system. We choose in each case the smallest 
coarse-graining volume for which transverse fluctuations capture 60% of the total variance of the fluctuations, i.e. 
||i5C"^|P > 0.6 ||<5C||^. This corresponds to Nb,, = 125 for the aging system and Nb^, = 216 for the equilibrium 
system. Table U contains a summary of the fitting parameters for those coarse graining volumes. 

Since the number of snapshots M used in the analysis cannot be increased indefinitely, it is important to test that 
the results are robust with respect to changes in M, and that meaningful results can be obtained even for relatively 
small values of M. There are two aspects to this question: one is to show that the phases 4>f{t) are well defined, by 
showing that their determination is robust with respect to changes in the number of snapshots M, and the other is 
to show that the magnitudes of the transverse, longitudinal and total fiuctuations are not singular as a function of 
M. With regards to the first aspect, we find that as the number of snapshots is increased, (^^{t) quickly converges 
to a relatively smooth function. Fig. ([3]) shows that (f>f{t) changes very little if the number of snapshots M used to 
compute it is 7 or larger in the case of the aging system or if it is 24 or larger in the case of the equilibrium system. It 
is noticeable that the minimum number Al^t of snapshots needed for a stable determination of 4>p{t) is comparable to 
the total variation A0 of the global phase over the interval considered: in the aging case Mgt ~ 7 and Acj) « 4, while 
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FIG. 2: Dependence of fluctuation strengths on the average number of particles Nb^ per coarse graining box, in the aging 
system with M — 100 and Nr = 100 (top panels) and in the system in equilibrium with M = 94 and A^^ = 100 (bottom panels). 
For each system, the left panel shows the magnitude \\dC\\ of the total fluctuations, the magnitude ||5C"^|| of the transverse 
fluctuations, and the magnitude ||5C^|| of the longitudinal fluctuations; and the right panels show the ratios ||(5C"^||/||5C|| and 
|i<5C^||/||5C||. The ratio ||5C^||/||(5C| increases and the ratio ||5C^||/||(5C|| decreases with increasing coarse graining size. For 
both systems, the first ratio is smaller than the second one at small coarse grainings, and larger than the second one at large 
coarse grainings. The two ratios cross at Nb,, ~ 30 for the aging system and at Nsp ~ 100 for the equilibrium system. 

in the equilibrium case Mgt ~ 24 and Acj) « 18. It is tempting to speculate that at least one or two configurations 
per relaxation time are needed in order to capture a minimum of information about the a relaxation in the system, 
and this would lead to the prediction that Mst/Acf) ^ 1-2. However at this point the data we have available are 
insufficient to draw any definite conclusions. 

We now consider the issue of the smoothness of the dependence of the magnitudes of the different kinds of fluctua- 
tions on the number M of snapshots. In Fig. (U) we plot the magnitude \ \SC\ \ of the total fluctuations, the magnitude 
[|(5C"^1| of the transverse fluctuations, and the magnitude |15C^|| of the longitudinal fluctuations as functions of 1/M. 
We find that all three fiuctuation magnitudes are gently increasing functions of M, and that each of them appears to 
approach a constant as M — > oo. For the specific coarse graining sizes shown in the figure, the longitudinal fiuctuations 
are always weaker than the transverse ones, but the opposite case is found for small enough coarse graining sizes. 

B. Local Phases and Local Relaxation Times 

Up to this point we have presented evidence in favor of the statement that the method we are introducing leads 
to a robust determination of the transverse and longitudinal components SC^ and SC^ of the fluctuations and of 
the fluctuating phases (j)r{t). We now use the method to show some examples of the phases (jipit) and their time 
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FIG. 4: Dependence of the magnitude ||<5C|| of the total fluctuations, the magnitude ||5C"^|| of the transverse fluctuations, and 
the magnitude ||5C^|| of the longitudinal fluctuations on the number A4 of snapshots used in their evaluation, for the aging 
system with Nb^ = 125 and Nr — 100 (left panel), and for the system in equilibrium with Nb,, ~ 216 and Nr = 100 (right 
panel). In all cases, the magnitude of each kind of fluctuation is a gently increasing function of AI that appears to extrapolates 
to a well-defined value as M — > oo. 



derivatives (l)p(t) = d4)p{t)/dt. The time derivatives are particularly interesting, not only because they are gauge 
invariant, but also because, according to Eq. (jlOp . they can in principle be interpreted as local relaxation rates or 
inverse local relaxation times Tp^{t). 

In Fig. ([S]) we show the time evolution of the local phases (j)f{t) at a fixed point r in the sample for different thermal 
histories, for both the aging system and the equilibrium system. For comparison, we also show the evolution of the 
global phase (j){t) , and we indeed find that for different thermal histories the local phase (j>f{t) gently fluctuates around 
the global value. 

In Figs. ^ and ([7]) we show the time evolution of the time derivative <j)f{t) of the local phase at the same point in 
space and for the same thermal histories as in the previous figure. In Fig. ([6]), (j>r{t) is rescaled so that the instantaneous 
global relaxation rate is unity, and results are shown both for the aging and the equilibrium cases. In the case of 
the equilibrium system, the rescaling is done approximately, by multiplying the time derivative of the phase by the 
a-relaxation rate Tq = l/dip/dt, where denotes a time average, and the fluctuating values of 4)p(t) are plotted as 
functions of the rescaled time t/xa, which is approximately equal to the global phase plus a constant. In the case of 
the aging system, the rescaling is done by dividing by d(j)/dt, and the scaled values of <j)f{t) are plotted as functions 
of the global phase 4i{t). In both cases the time derivative of the local phase fluctuates strongly, often becoming more 
than twice as large as the global relaxation rate d(j)/dt, or even becoming negative. It should be pointed out that 
at the times when 4'f{t) becomes negative, its interpretation as a local relaxation rate 1/Tf{t) becomes problematic. 
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FIG. 5: Comparison between the local phase <j>r{t) as a function of time at a fixed point in space for different thermal histories 
and the global phase (j}{t). Shown for the aging system, with Nsf, = 125 and M = 100 (left panel), and for the equilibrium 
system, with A''^^ — 216 and A/ = 94 (right panel). 



However, from a more general point of view, it is not surprising that transient fluctuations in a small region of the 
sample may give rise to changes that make the configuration of the region temporarily become closer to what it was in 
the past, thus making the two time correlation Cp(t,t') increase over time instead of decreasing, with the effect that, 
for some time at least, (pp{t) becomes negative. In Fig. ([7]) we plot the time derivative of the local phase in the aging 
system, without rescaling it, as a function of time in Lennard- Jones units. For times longer than 2.5 Lennard- Jones 
time units, the global relaxation rate d<j)/dt decreases as a function of time, which is a direct consequence of the aging 
in the system, namely the fact that the dynamics becomes slower as the relaxation progresses. In this figure we can 
notice again the strong fluctuations of the time derivative of the local phase with respect to the global value, but we 
also notice that both the typical value and the fluctuations of (f)p{t) decrease together with the decrease of d(j>/dt. Also 

shown in the same figure is a fit of the global relaxation rate by the form dcjj/dt — , where to 3.4 x 10^'* and 

a ~ —0.91, a result that is in agreement with the results in Refs. j26l. [27j. 




FIG. 6: Time evolution of the time derivative of the local phases (pp, rescaled so that the time derivative of the global phase 
is unity. In the left panel, corresponding to the aging system with A''^^ — 125 and M — 100, the time derivative of the local 
phases is rescaled by dividing by the global time dependent relaxation rate and it is plotted as a function of the global 
phase (;/>(t). In the right panel, corresponding to the equilibrium system with Nsp ~ 216 and M = 94, the rescaling of (pp is 
performed approxmately by multiplying it with the a-relaxation rate = l/d(f>/dt, and the rescaled time derivative is plotted 
as a function of the rescaled time t/ra, which is approximately equal to the global phase (pit) plus a constant. 
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FIG. 7: Time evoluion of the time derivative of the local phase, (j)^, for the aging system with A''^^ — 125 and M = 100, with 
no reseating. The global relaxation rate is also shown, together with a fit of the form d(j)/dt = (t/to)", where to ~ 3.4 x 10"'' 
and a ^ —0.91. 



V. SUMMARY AND CONCLUSIONS 



In this work we have shown that a quantitative description of glassy relaxation in terms of local fluctuating phases 
and relaxation times is possible. On one hand, starting from the case of exponential relaxation we presented a line 
of phenomenological arguments that shows how a local phase function can emerge in the description of the 

data, and how its time derivative can, under certain circumstances, be interpreted as a local fluctuating relaxation 
rate 1/Ti^{t). On the other hand, we showed that in a theoretical framework that postulates the presence of a broken 
symmetry under time reparametrizations, it is expected that as the dynamics becomes more glassy, fluctuations should 
be dominated by Goldstone modes - i. e. transverse fluctuations - which are described naturally in terms of the same 
local fluctuating phases (f>p{t). Besides establishing a connection between the two points of view, we have also derived 
a practical method for extracting the phases (j^pit) from numerical or experimental data. 

We have applied this method to two numerical simulation datasets, one corresponding to a system undergoing 
aging at a temperature slightly below the mode coupling critical temperature and another one corresponding to a 
system in equilibrium at a temperature slightly above T^. We have found that in both cases the results are robust 
with respect to changes in the number M of position snapshots used in the analysis, as long as M is not too small. 

The time reparametrization symmetry framework predicts that for larger coarse graining regions and lower tem- 
peratures the transverse fluctuations should be dominant, and indeed, this is what we found in our results. For both 
systems, the ratio between the magnitude of the transverse fluctuations and the magnitude of the total fluctuations 
grows monotonously with the coarse graining size. For the smallest coarse graining sizes, longitudinal fluctuations 
dominate, but the transverse fluctuations dominate for large coarse graining sizes. Also, as expected, the range of 
coarse graining sizes for which transverse fluctuations dominate is wider in the system that is in contact with a heat 
reservoir at a lower temperature than in the system that is in equilibrium at a higher temperature. 

We have shown more detailed results for the fluctuating phase (j)f{t) and its time derivative <j)p{t), for a coarse 
graining size such that the transverse fluctuations capture slightly more than 60% of the total fluctuations. For this 
coarse graining size, we found that both in the equilibrium system and in the aging system, the time derivative of 
the local phase fluctuates rather strongly over timescales of the order of the instantaneous global relaxation times 
or shorter. These fluctuations are strong enough that the time derivative often takes negative values. The presence 
of those negative values, although not surprising, makes an interpretation of 4>r{t) as a relaxation rate somewhat 
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problematic. It is tempting to speculate that this interpretation may be more cleanly applicable for large enough 
coarse graining sizes, for which the total fluctuations are expected to be weak enough that (t>p{t) should always remain 
positive. 

Since our method involves fitting the local two-time correlations Cp{t,t') by the expression g[(j)p(t) — (j)p(t')], one 
of our results is the functional form for g{x). In our fits we modeled this function as a stretched exponential g{x) — 
g_EA exp(— Phenomenological arguments have been made that would indicate that the relaxation should be 
exponential for small enough coarse graining regions. If this was the case, we should find the fitting parameter (3 
approaching unity for smaller coarse graining sizes. Our results in this respect are inconclusive: although the fitting 
parameter /3 does approach unity as the coarse graining size is reduced, and even becomes larger than unity in the case 
of the aging system, this happens for a range of coarse graining sizes in which the fluctuations are mostly longitudinal, 
or in other words most of the fluctuations are "explained" by the fitting residuals and not by the fit itself. 

One of the ways in which the study of glasses is more difficult than the study of other condensed matter systems 
is that in the case of glasses there are no simple observables that allow probing glassiness directly. For example, 
if we want to probe magnetic ordering, a one-time observable, the magnetization (S{f,t)), is enough to detect the 
presence of ordering, and magnetic fluctuations can be probed by the two-point function {S{f,t) ■ S(r',t')). For 
the case of crystalline ordering, the ordering can be probed by the one-time Fourier transformed particle density 
p{k,t) = iV^^ ^^■(exp[iA: ■ rj{t)]), and correlations can be probed by the dynamic structure factor S{k,t — t') = 

N^^ X]j i' (s^pI^*^ ■ 'S'lO] exp[ifc • rji{t')]). In the case of structural glasses, by contrast, the presence of a glass state 
is normally detected by probing the freezing of the density fluctuations, which already requires a two-time observable 
such as the intermediate scattering function Fs{q^ t, t') or the mean square displacement A(t, t'); and the corresponding 
fluctuations are quantified by a 4-point correlation function such as Gi{r^t — t') or its Fourier transform S4{q,t,t'), 
which correspond to taking products of two-time observables probing different spatial points in the system. The 
description of glassy dynamics that we have presented here, although it is extracted from the two-time quantities that 
probe freezing in the system, has the advantage that it is expressed in terms of one time quantities, the local phases 
4>r{i), or, more properly, their gauge invariant time derivatives 4>fit)- This should make it possible, in principle, to 
probe correlations in the dynamics by computing two-point functions, such as {(j)p{t)(j)n{t')), as opposed to four-point 
functions. One way of thinking about this is that this advantage in simplicity is gained by coarse graining both in 
space and in time. The coarse graining in space is explicit, and has been discussed extensively in this work; the coarse 
graining in time happens in the sense that the fitting procedure at the heart of the method focuses on extracting 
information about the slow part of the dynamics and ignores the fast part. On one hand, it could be argued that this 
coarse graining comes with two kinds of costs: one is the loss of some information about processes that are "fast" in 
space and time, and the other is the fact that some extra steps have been introduced in the analysis of the data. On 
the other hand, it is possible that this loss of information and these extra steps in the analysis may make it easier to 
observe more cleanly the essential features of the dynamics. 

At this point we have introduced a new general analysis tool to study glassy dynamics. Since the starting point 
of this analysis is the determination of local two-time correlations, it can be applied to any experiment or numerical 
simulation data set which provides enough information to compute those correlations. This includes not only results 
from numerical simulations, but also results from visualization experiments in which the particle positions are 

determined as a function of time. More generally, there are other kinds of experiments, such as photon correlation 
imaging [sl], which provide data for local two-time correlations, and in this case too this type of analysis is applicable. 

The investigation presented here leaves many questions open to be addressed. Those questions include the deter- 
mination of the basic statistical properties of the fluctuating quantity (j)p{t) , such as its probability distribution and 
its correlation functions, the study of how those properties vary with changes in the control parameter, namely either 
the temperature or the volume fraction, and also the study of how they change across different systems exhibiting 
glass-like dynamics. In particular, the correlation function {(j)p{t)(j) ~,{t')) encodes important information about dy- 
namical heterogeneity: in the equal-time t — t' case it gives very direct information about the spatial extent of the 
heterogeneous regions, and in the equal-point r — r' case it gives information about the persistence in time of those 
regions. Another set of questions that would be worth investigating is the relationship between the present descrip- 
tion of the dynamics, and descriptions based on other probes such as the dynamical susceptibility Xi or the various 
point-to-set correlation functions. Another set of questions that could be addressed is a more detailed study of the 
form of the function g{x), and in particular the question of the possible presence of simple exponential relaxation at 
relatively small coarse graining scales. 
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Appendix A: Derivation of the Method 



Here we present some extra details about the derivation of the method for extracting local phases from numerical 
or experimental data, which we omitted in Sec. IIIII We start by considering the minimization of e{{4>i, • ■ • , </>m}; cS), 
defined by Eq. (1^51) . By using the definition of in Eq. (P7|) . and explicitly writing the terms in the Taylor expansion 
g''^^{4'i ~ 4'u o^) we rewrite e to obtain 

• • • , 0m}; d) = r\M) (ff ('J'r- ") - {9 ('i'j.; a) + g' (<i>r.; «) [{<f>j - -^z) - + •■•})' 

l<i<j<M 

= E {9'i'^>,^■,d)[icb,-cb.)-^,.] + ■■■f (Al) 

' l<i.j<M 

where the terms omitted are of quadratic order in — 4>j) ~ )• In ttis last line we have allowed both i > j and 
i < j terms to appear in the sum, and we have defined = so that i = j terms are identically zero. We now 
impose the condition that all derivatives of e with respect to (jfk be zero, 

de 

where we have used the fact that g{—x) = g{x) and therefore g'{—x) — —g'{x). Hence for 1 < fc < M, we have the 
minimization conditions: 

M M M 

i—1 i—1 i—1 

We express the minimization conditions in terms of the equivalent matrix equation 

li) = A(j3, (A4) 

with the definition of the M-component vector w 

M 

t&fc = ^.g'^ (A5) 

and the M x M matrix A 

M 

Au = 5" {<^^k ; d) - 4» 5" ; ■ 

The system of equations described by Eq. (|A4[) is singular, since the nonzero vector (^o = (I7 1) satisfies A(t)Q — 0, 
which indicates that the rank of the matrix A is at most M — 1. In particular, any solution of Eq. (|A4I) can be 
modified by adding an arbitrary multiple of to obtain another solution. This degeneracy is associated with the 
gauge symmetry that the problem has under the shift (pk 4'k+ P, which is a consequence of the fact that in Eq. (fT2)) 
the phases appear in a difference. To eliminate this degeneracy, we impose the additional condition that the projection 
of the solution (j) on the direction of (f>Q should be zero, thus obtaining the gauge fixing condition of Eq. (1281) . By 
combining the gauge fixing condition with Eq. (jA4p above, we obtain a system of Af + 1 equations with M unknowns: 

w ^ A(f and 00 • = 0. (A7) 

Since one of the equations is a linear combination of the others, the solution is well defined and unique. Hence, we 
omit the first equation, for fc = 1, which is redundant, and solve 



(A8) 



where w and A are given by 



_ ^ I : fc = 1 

~ \wk^ YZi ff" i'^^k-, 5) : fc ^ 1 

A = \ ^ : fc = 1 , . „x 

" \Au= .g'2 ; d) - 6k. ^ -ii 9" i'^ik ; cf) : fc ^ 1. ^ ' 
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which is the same as Eq. ((3T|) . 

We now turn to the determination of the local phase fluctuations Scj^ip, which requires minimizing the expression 
for e in Eq. ([22|) . Here the derivation is very similar to the previous case. We start by using the definition of SCp{t, t') 
in Eq. ([T5|) to rewrite e as 

e{{5(^,^{U)}-a,{C,{tj,U)}) = ^ {5Cr-{tj,U) - g' m,) - cj,{U);d) [SMW " SMU)]f ■ (AlO) 



We then impose the condition that the derivatives with respect to all of the phase fluctuations S(f>kr must be zero, 

= — ^ = -2t^-\M) V [SCk^f - a' {C^k - (t'^■,d) {S<f>kr " S^,p)] q' {<f>k ' d) . (All) 
oScpkr 

This is equivalent to imposing, for 1 < fc < M, the conditions 

M M M 

g' {<l>k - <l>i;d) SCkif = Hkr^^g''^ {<l>k - <l>i-d) g"^ {(j>k - 4ii;d) 5<j)ip. (A12) 

i—1 i—1 i—1 

This set of conditions can be written as the matrix equation 

wp = A5(j)p, (A13) 

where 

M M 

Wkr='^g'{4'k-'4'j;d)SCkjf and Aki = Skt^g''^ {(l)k ~ <l>j;d) - g''^ {(j>k ~ (j>i;d) . (A14) 

From here, the steps are almost identical to those for the determination of the global phases. Just as in that case, we 
have a degeneracy resulting from the shift symmetry under the transformation Scji^p 5(j)kp + pp^ where pp is time 
independent. We again modify the system of equations by imposing an additional constraint, given by Eq. p2p . that 
removes the degeneracy. Here again we omit the redundant equation corresponding to fc = 1, and thus finally recover 
Eqs. (P I) and ([^S ]) . 

Appendix B: Proof of the Orthogonality of 5C^ and 5C^ 

In this section we prove Eq. (HDl) . We start by applying the definition of the Euclidean scalar product in Eq. 

(5C^|<5C^) = -— |-- T.(^CTit,,U)SCi{t„U)). (Bl) 

^ ' ' l<i<i<M r 

We now insert the definitions of 5C^ and 5G^ , given by Eqs. (fT7|) and ([T9l) . and obtain 

^ ^ ' l<j,i<M r 

where we have removed the restriction j > i. By using the symmetry of the terms under the exchange of indices 
i ■(^ j, we rewrite this expression as 

/ M M \ 

{5C^\5C^) = ^^^^^^ _ 50,^.^ g'(0, - </..) {5C,,p ~ g'{<i,, - 0,; a) {8<i>,p - Jc/..^}^ • (B3) 

We now recognize that, by Eq. (jAllI) . the sum over the index i is zero, and therefore the whole expression is zero, 
which gives us the result that (i5C-^|(5C^) = 0. 



[1] P. G. Debenedetti and F. H. Stillinger, Nature 410, 259 (2001). 



16 



[2] M. D. Ediger, Annu. Rev. Phys. Chem. 51, 99 (2000). 
[3] H. Sillescu, J. Non-Crystal. Solids 243, 81 (1999). 

[4] E. Vidal Russell, N. E. Israeloff, L. E. Walther, and H. Alvarez Gomariz, Phys. Rev. Lett. 81 1461 (1998); L. E. Walther, 
N. E. Israeloff, E. Vidal Russell, and H. Alvarez Gomariz Phys. Rev. B 57 R15112 (1998); E. Vidal RusseU and N. E. 
Israeloff, Nature 408, 695 (2000). 

[5] E. R. Weeks, J. C. Crocker, A. C. Levitt, A. B. Schofield, and D. A. Weitz, Science 287, 627 (2000); 

[6] E. R. Weeks, and D. A. Weitz, Phys. Rev. Lett. 89 095704 (2002). 

[7] R. E. Courtland and E. R. Weeks, 2003 J. Phys.: Condens. Matter 15 S359. 

[8] L. Cipelletti, H. Bissig, V. Trappe, P. Ballesta, and S. Mazoyer, J. Phys.: Condens. Matter 15, S257 (2003). 

[9] A. S. Keys, A. R. Abate, S. C. Glotzer, and D. J. Durian, Nature Physics 2007, 260-264, (2007). 
[10] W. Kob, C. Donati, S. J. Phmpton, P. H. Poole, and S. C. Glotzer, Phys. Rev. Lett. 79, 2827 (1997). 
[11] N. Lacevic, F. W. Starr, T. B. Schroder and S. C. Glotzer, J. Chem. Phys. 119, 7372 (2003). 
[12] C. TonineUi, M. Wyart, L. Berthier, G. Biroh, and J.-P. Bouchaud, Phys. Rev. E 71, 041505 (2005). 
[13] J. P. Garrahan and D. Chandler, Phys. Rev. Lett. 89, 035704 (2002). 
[14] D. Chandler, and J. P. Garrahan, Annu. Rev. Phys. Chem. 61, 191 (2010). 
[15] X. Xia, P. G. Wolynes, Phys. Rev. Lett. 86, 5526 (2001). 
[16] V. Lubchenko and P. G. Wolynes, Annu. Rev. Phys. Chem. 58, 235 (2007). 

[17] C. Chamon, M. P. Kennett, H. E. Castillo, and L. F. Cughandolo, Phys. Rev. Lett. 89, 217201 (2002). 

[18] H. E. Castillo, C. Chamon, L. F. Cughandolo, and M. P. Kennett, Phys. Rev. Lett. 88, 237201 (2002). 

[19] H. E. Castillo, C. Chamon, L. F. Cughandolo, J. L. Iguain, and M. P. Kennett, Phys. Rev. B 68, 134442 (2003). 

[20] H. E. Castillo and A. Parsaeian, Nature Physics 3, 26 (2007). 

[21] H. E. Castillo, Phys. Rev. B 78, 214430 (2008); 

[22] A. Parsaeian and H. E. Castillo, Phys. Rev. E 78, 060105(R) (2008). 

[23] A. Parsaeian and H. E. Castillo, arXiv:0811.3190 (2008). 

[24] A. Parsaeian and H. E. Castillo, Phys. Rev. Lett. 102, 055704 (2009). 

[25] G. A. Mavimbela, H. E. Castillo, J. Stat. Mech. (2011) P05017. 

[26] K. E. Avila, H. E. Castillo, and A. Parsaeian, Phys. Rev. Lett. 107, 265702 (2011). 

[27] K. E. Avila, H. E. Castillo, and A. Parsaeian, in preparation. 

[28] K. Binder, W. Kob, Glassy Materials and Disordered Solids (World Scientific Publishing Co. Pte. Ltd, 2005). 
[29] R. Bohmer et. al, J. of Non-Crystalhne Solids, 235 (1998). 
[30] L. F. Cughandolo, J. Kurchan, 1994 J. Phys. A 27 5749 

[31] L. F. Cughandolo, J. Kurchan, 1993 Phys. Rev.Lett. 71(1) 173 

[32] S. Franz, G. Parisi, F. Ricci-Tersenghi and T. Rizzo, arXiv:100 fl74 6l (2010); larXiv:1008.0996l (2010). 
[33] S. Franz, G. Parisi, F. Ricci-Tersenghi and T. Rizzo, Eur. Phys. J. E 34, 102 (2011). 
[34] W. Zwerger, Phys. Rev. Lett. 92, 027203 (2004). 

[35] A. Duri, D. A. Sessoms, V. Trappe, and L. Cipelletti, Phys. Rev. Lett. 102, 085702 (2009) 



